RGI bwt output
library(readxl)
setwd("/Users/karlavasco/Library/CloudStorage/OneDrive-MichiganStateUniversity/Manning_lab/Mastitis_project/metagenome_KV/tables/resistome/")
resistome_full <- read_excel("rgi_mastitis.xlsx", sheet = "RGI_bwt_allele_seqname_noheader")
Filter mutations and eliminate them from the final resistome table
library(data.table)
library(dplyr)
`%notin%` <- Negate(`%in%`)
mutants1 <- resistome_full %>% filter(`Reference Sequence` %like% "conferring")
mutants2 <- resistome_full %>% filter(`Reference Sequence` %like% "muta")
resistome <- resistome_full %>% filter(`Reference Sequence` %notin% c(mutants1$`Reference Sequence`, mutants2$`Reference Sequence`))
Importing sample data
setwd("/Users/karlavasco/Library/CloudStorage/OneDrive-MichiganStateUniversity/Manning_lab/Mastitis_project/metagenome_KV/tables/metadata/")
data <- read_excel("IMM_ceftiofur_metadata.xlsx", sheet = "metadata")
metagenome_stats <- read_excel("IMM_ceftiofur_metadata.xlsx", sheet = "metagenome_stats")
temperature <- read_excel("IMM_ceftiofur_metadata.xlsx", sheet = "temperature")
diet <- read_excel("IMM_ceftiofur_metadata.xlsx", sheet = "diet")
Merging tables as metadata
metadata <- merge(data, diet) %>%
merge(temperature) %>%
merge(metagenome_stats) %>%
merge(resistome) %>%
mutate(Treatment = factor(Treatment, levels = c("Control", "Antibiotic")),
Time_tx = factor(Time_tx, levels = c("Day -1", "Week 1", "Week 5", "Week 9")))
New column with normalized abundance based on Genome equivalents
metadata <- metadata %>% mutate(normalized_abundance = Depth/genome_equivalents*100)
Making table
abundance_feature <- metadata %>%
select(normalized_abundance, Treatment, Time_tx, sample_ID) %>%
group_by(Treatment, Time_tx, sample_ID) %>%
summarise(Abundance = sum(normalized_abundance)) %>% as.data.frame() %>% mutate(Type = "Resistome")
abundance_feature$Treatment <- factor(abundance_feature$Treatment, levels= c("Control","Antibiotic"))
Abundance plot
library(rstatix)
library(ggpubr)
#Calculating p-values between treatments by time point
stat.test <- abundance_feature %>%
filter(sample_ID %notin% c("MA028.7")) %>%
group_by(Time_tx) %>%
wilcox_test(Abundance ~ Treatment, alternative = "less", paired = T) %>%
adjust_pvalue(method = "bonferroni") %>%
add_significance("p.adj") %>%
add_xy_position(x = "Time_tx", dodge = 0.8)
stat.test
Abundance_boxplot <- abundance_feature %>%
ggboxplot(x = "Time_tx", y = "Abundance",
color = "Treatment", palette = "jama", fill = "Treatment", alpha = 0.5,
add = c("jitter"), outliers = F,notch = F, outlier.shape = NA,
facet.by = "Type",
xlab = F, ylab = "Abundance score") +
stat_pvalue_manual(stat.test, label = "p", tip.length = 0,
y.position = c(81.55,86.24,71.309,75)) +
ylim(30,90)
Abundance_boxplot
Table with abundance
abundance_class <- metadata %>%
dplyr::select(normalized_abundance, `Drug Class`, Time_tx, Treatment) %>%
dplyr::group_by(`Drug Class`) %>%
dplyr::summarise(Proportion = sum(normalized_abundance*100/sum(metadata$normalized_abundance))) %>%
as.data.frame() %>%
mutate(Type = "Resistome") %>%
arrange(desc(Proportion))
top_20 <- abundance_class %>% dplyr::select(`Drug Class`,Proportion) %>% top_n(10)
abundance_class$Proportion <- format(round(abundance_class$Proportion, 2), nsmall = 2)
abundance_class$Proportion <- as.numeric(as.character(abundance_class$Proportion))
head(abundance_class)
Drug-class plot
library(stringr)
library(ggsci)
class_plot <- abundance_class %>%
filter(`Drug Class` %in% top_20$`Drug Class`) %>%
mutate(across('Drug Class', str_replace, ' antibiotic', '')) %>%
mutate(across('Drug Class', str_replace, ' antibiotic;', ';')) %>%
ggplot(aes(y = reorder(`Drug Class`, +Proportion),
x = Proportion,
color =`Drug Class`,fill=`Drug Class`,
ylab=F,
label = Proportion)) +
geom_bar(stat="identity") +
theme_minimal() +
geom_text(size = 3, color = "black", hjust = -.2, position = position_identity())+
theme(legend.position = "none", axis.title.y = element_blank()) +
scale_fill_simpsons() +
scale_color_simpsons() +
xlim(0,48.5)
class_plot
Saving plot
setwd("/Users/karlavasco/Library/CloudStorage/OneDrive-MichiganStateUniversity/Manning_lab/Mastitis_project/metagenome_KV/figures")
ggsave(plot=class_plot, "resistome_class.png", width = 6.9, height = 3.2,bg="white")
Calculating the number of alleles to manually add in Drug-class plot
number_alleles_class <- metadata %>%
dplyr::select(`ARO Term`, `Drug Class`) %>%
distinct() %>%
mutate(count = 1) %>%
dplyr::group_by(`Drug Class`) %>%
summarise(total = sum(count)) %>%
arrange(desc(total))
number_alleles_class
Selecting the top 70 args
`%notin%` <- Negate(`%in%`)
presence <- metadata %>%
select(`ARO Term`, sample_ID) %>%
group_by(`ARO Term`) %>%
tally() %>% arrange(desc(n))%>%
filter(n>=15)
presence
abundance_allele <- metadata %>%
dplyr::select(normalized_abundance, `ARO Term`, `Reference Model Type`) %>%
dplyr::rename(Allele = `ARO Term`) %>%
dplyr::filter(Allele %in% presence$`ARO Term`) %>%
dplyr::group_by(Allele) %>%
dplyr::summarise(Abundance = sum(normalized_abundance)) %>% as.data.frame()
top70 <- top_n(abundance_allele, 70) %>% arrange(desc(Abundance))
Matrix with sequence_id as columns and features as rows
library(tibble)
top70_mx <- metadata %>% select(sample_ID,`ARO Term`,normalized_abundance)%>%
dplyr::select(sample_ID,`ARO Term`,normalized_abundance) %>%
dplyr::rename(Allele = `ARO Term`) %>%
dplyr::filter(Allele %in% top70$Allele) %>%
dplyr::group_by(sample_ID, Allele) %>%
dplyr::summarise(Abundance = sum(normalized_abundance)) %>%
tidyr::spread(Allele, Abundance) %>%
merge(metadata, by="sample_ID") %>%
dplyr::arrange(Time_tx,Treatment) %>%
dplyr::select(sample_ID, top70$Allele) %>%
distinct() %>%
remove_rownames %>% column_to_rownames(var="sample_ID") %>%
as.matrix(rownames=TRUE)
top70_mx[is.na(top70_mx)] <- 0
Annotations for plot
annotations <- metadata %>%
select(sample_ID,Treatment, Time_tx) %>%
distinct() %>%
remove_rownames %>% column_to_rownames(var="sample_ID")
anno_color <- list(Treatment = c(Control = "#374E55FF",
Antibiotic = "#DF8F44FF"),
Time_tx = c(`Day -1` = "#3B4992FF",
`Week 1` = "#EE0000FF",
`Week 5` = "#008B45FF",
`Week 9` = "#631879FF"))
Heatmap plot
library(pheatmap)
library(viridis)
# Gene names in italics
newnames <- lapply(
colnames(top70_mx),
function(x) bquote(italic(.(x))))
# Heatmap
heatmap <- pheatmap(
mat = t(log10(top70_mx+0.0000001)),
border_color = NA,
show_colnames = F,
show_rownames = T,
angle_col = 90,
color = viridis(100), annotation_col = annotations,
annotation_colors = anno_color,
annotation_names_row = F,
cluster_cols = F,
cluster_rows = T,
clustering_method = "ward.D",
gaps_row = FALSE,
labels_row = as.expression(newnames)
)
heatmap
Saving heatmap
setwd("/Users/karlavasco/Library/CloudStorage/OneDrive-MichiganStateUniversity/Manning_lab/Mastitis_project/metagenome_KV/figures")
ggsave(plot=heatmap, "heatmap_resistome.png", width = 9, height = 10)
Feature table used for alpha diversity: Matrix with sequence_id as columns and features as rows
library(tibble)
feature_alpha <- metadata %>% select(sequence_id,`Reference Sequence`,Depth) %>%
spread(sequence_id, Depth) %>%
remove_rownames %>% column_to_rownames(var="Reference Sequence") %>%
as.matrix(rownames=TRUE)
feature_alpha[is.na(feature_alpha)] <- 0
Taxonomy table
taxonomy <- resistome_full %>% select(`Reference Sequence`, `Reference Model Type`, `Drug Class`, `Resistance Mechanism`,`AMR Gene Family`,`ARO Term`) %>% distinct() %>%
rename(Model = `Reference Model Type`,Class=`Drug Class`,Mechanism=`Resistance Mechanism`,Family=`AMR Gene Family`,Allele=`ARO Term`) %>%
remove_rownames %>% column_to_rownames(var="Reference Sequence") %>%
as.matrix()
Metadata to matrix
metadata_matrix <- merge(data,metagenome_stats) %>% remove_rownames %>% column_to_rownames(var="sequence_id")
Phyloseq object
library("phyloseq")
feature_alpha <- feature_alpha*10^4
mode(feature_alpha) <- "integer"
feature_alpha[is.na(feature_alpha)] <- 0
#import as phyloseq objects
OTU = otu_table(feature_alpha,taxa_are_rows=TRUE)
TAX = tax_table(taxonomy)
META = sample_data(metadata_matrix)
physeq=phyloseq(OTU,TAX,META)
alpha_diversity <- estimate_richness(physeq, measures = c("Shannon", "Observed"))
df_alpha <- data.frame(alpha_diversity, sample_data(physeq))
alpha_table <- reshape2::melt(df_alpha, measure.var=c("Shannon","Observed"),id.vars=c("Time_tx","Treatment","sample_ID","study_ID","Block","Cohort")) %>%
rename(Index = variable)
alpha_table$value = as.numeric(alpha_table$value)
alpha_table$Treatment <- factor(alpha_table$Treatment, levels= c("Control","Antibiotic"))
alpha_table <- alpha_table %>% mutate(Type = "Resistome")
Shannon diversity
#Calculating p-values between treatments by time point
stat.test <- alpha_table %>%
filter(Index == "Shannon", sample_ID %notin% c("MA028.7")) %>%
group_by(Time_tx) %>%
wilcox_test(value ~ Treatment, alternative = "greater", paired = T) %>%
adjust_pvalue(method = "bonferroni") %>%
add_significance("p.adj") %>%
add_xy_position(x = "Time_tx", dodge = 0.8)
stat.test
#Boxplot of Shannon diversity between treatment groups over time
shannon_boxplot = alpha_table %>%
filter(Index == "Shannon") %>%
ggboxplot(x = "Time_tx", y = "value",
color = "Treatment", palette = "jama", fill = "Treatment", alpha = 0.5,
add = c("jitter"), notch = F, outlier.shape = NA, facet.by = "Type",
xlab = F, ylab = "Shannon Index", legend = "none") +
ylim(3.1,5.1)+
stat_pvalue_manual(stat.test, label = "p", tip.length = 0)
shannon_boxplot
Observed diversity
#Calculating p-values between treatments by time point
stat.test <- alpha_table %>%
filter(Index == "Observed", sample_ID %notin% c("MA028.7")) %>%
group_by(Time_tx) %>%
wilcox_test(value ~ Treatment, alternative = "greater", paired = T) %>%
adjust_pvalue(method = "bonferroni") %>%
add_significance("p.adj") %>%
add_xy_position(x = "Time_tx", dodge = 0.8)
stat.test
library(ggpubr)
library(ggsci)
#Boxplot of Observed diversity between treatment groups over time
Observed_boxplot = alpha_table %>%
filter(Index == "Observed") %>%
ggboxplot(x = "Time_tx", y = "value",
color = "Treatment", palette = "jama", fill = "Treatment", alpha = 0.5,
add = c("jitter"), notch = F, outlier.shape = NA, facet.by = "Type",
xlab = F, ylab = "Shannon Index", legend = "none") +
ylim(90,400)+
stat_pvalue_manual(stat.test, label = "p", tip.length = 0)
Observed_boxplot
Feature table used for beta diversity and differential abudance analyses: Matrix with sequence_id as columns and features as rows
library(tibble)
feature_table <- metadata %>% select(sequence_id,`Reference Sequence`,normalized_abundance) %>%
spread(sequence_id, normalized_abundance) %>%
remove_rownames %>% column_to_rownames(var="Reference Sequence") %>%
as.matrix(rownames=TRUE)
feature_table[is.na(feature_table)] <- 0
Phyloseq object for beta diversity
library("phyloseq")
#import as phyloseq objects
OTU = otu_table(feature_table,taxa_are_rows=TRUE)
#(tree was already imported as a phyloseq object)
physeq=phyloseq(OTU,TAX,META)
Agglomerating ARGs by allele
physeq <- tax_glom(physeq, "Allele")
PERMANOVA by Time
library(microbial)
beta <-betatest(physeq,group="Time_tx")
beta
Bray-Curtis plot
bray_resistome <- plotbeta(
physeq,
group="Time_tx",
shape = "Treatment",
distance = "bray",
method = "PCoA",
color = F,
size = 3,
ellipse = F) +
labs(color = "Time", shape = "Treatment", fill="Time") +
annotate("text", x = 0, y = -0.22, label = expression(paste("PERMANOVA, ",F ,"= 11.98, ",paste(italic('p')),"=0.001")), colour = "black", size = 4) +
ggtitle("Resistome") +
stat_ellipse(geom = "polygon",
aes(fill = Time_tx, group=Time_tx),
alpha = 0.1)+
scale_fill_aaas() +
scale_color_aaas()+
theme_bw()
bray_resistome
Arranging microbiome plots
resistome_plots <- ggarrange(Abundance_boxplot, class_plot, shannon_boxplot, bray_resistome, labels = c("A","C","B","D"), widths = c(1,2,1,2), heights = c(1.2,1), nrow = 2,ncol = 2)
resistome_plots
Saving plot
setwd("/Users/karlavasco/Library/CloudStorage/OneDrive-MichiganStateUniversity/Manning_lab/Mastitis_project/metagenome_KV/figures")
ggsave(plot=resistome_plots, "Fig5_resistome_plots2.png", width = 10, height = 6, bg="white")
List of betalactamases with cephalosporin activity
library(data.table)
taxonomy_df <- as.data.frame(taxonomy)
betalactamases <- taxonomy_df[taxonomy_df$Family %like% "beta-lactamase", ]
esbls <- betalactamases[betalactamases$Class %like% "ceph", ]
ceph <- taxonomy_df[taxonomy_df$Class %like% "ceph", ]
ESBLs abundance matrix
betalactamases_mx <- metadata %>%
select(sample_ID,`ARO Term`,normalized_abundance) %>%
dplyr::select(sample_ID,`ARO Term`,normalized_abundance) %>%
dplyr::rename(Allele = `ARO Term`) %>%
dplyr::filter(Allele %in% esbls$Allele) %>%
dplyr::group_by(sample_ID, Allele) %>%
dplyr::summarise(Abundance = sum(normalized_abundance)) %>%
tidyr::spread(Allele, Abundance) %>%
merge(metadata, by="sample_ID") %>%
dplyr::arrange(Time_tx,Treatment) %>%
dplyr::select(sample_ID, esbls$Allele) %>%
distinct() %>%
remove_rownames %>% column_to_rownames(var="sample_ID")%>%
as.matrix(rownames=TRUE)
betalactamases_mx[is.na(betalactamases_mx)] <- 0
ESBLs Heatmap
newnames <- lapply(
colnames(betalactamases_mx),
function(x) bquote(italic(.(x))))
#HEATMAP RA USING ln color pallet RBrewer
heatmap <- pheatmap(
mat = t(log10(betalactamases_mx+0.0000001)),
border_color = NA,
show_colnames = F,
show_rownames = T,
angle_col = 90,
color = viridis(100),
annotation_col = annotations,
annotation_colors = anno_color,
annotation_names_row = F,
cluster_cols = F,
cluster_rows = F,
clustering_method = "ward.D",
gaps_row = FALSE,
labels_row = as.expression(newnames)
)
heatmap
Saving ESBL heatmap
setwd("/Users/karlavasco/Library/CloudStorage/OneDrive-MichiganStateUniversity/Manning_lab/Mastitis_project/metagenome_KV/figures")
ggsave(plot=heatmap, "heatmap_esbls.png", width = 9, height = 18)
Cephalosporin resistance abundance table
cephalosporin <- metadata %>%
dplyr::select(sample_ID,`ARO Term`,normalized_abundance,Time_tx,Treatment) %>%
dplyr::rename(Allele = `ARO Term`) %>%
dplyr::filter(Allele %in% ceph$Allele) %>%
dplyr::group_by(sample_ID, Time_tx,Treatment) %>%
dplyr::summarise(Abundance = sum(normalized_abundance)) %>%
dplyr::mutate(Type="Cephalosporin resistance")
cephalosporin$Treatment <- factor(cephalosporin$Treatment, levels= c("Control","Antibiotic"))
Cephalosporin resistance plot
#Calculating p-values between treatments by time point
stat.test <- cephalosporin %>%
filter(sample_ID %notin% c("MA028.7")) %>%
group_by(Time_tx) %>%
wilcox_test(Abundance ~ Treatment, alternative = "less", paired = T) %>%
adjust_pvalue(method = "bonferroni") %>%
add_significance("p.adj") %>%
add_xy_position(x = "Time_tx", dodge = 0.8)
stat.test
cephalosporin_boxplot = cephalosporin %>%
ggboxplot(x = "Time_tx", y = "Abundance",
color = "Treatment", palette = "jama", fill = "Treatment", alpha = 0.5,
add = c("jitter"), notch = F, outlier.shape = NA,
facet.by = "Type",
xlab = F, ylab = "Abundance score") +
ylim(2,31)+
stat_pvalue_manual(stat.test, label = "p", tip.length = 0) +
geom_signif(aes(y = Abundance, x = Time_tx),
comparisons = list(c("Day -1","Week 1"),c("Week 5", "Week 9")),
test = "wilcox.test",
test.args=list(alternative = "less", var.equal = FALSE, paired=FALSE),
color = "black", tip_length = 0.02, y_position = c(21,23)) +
geom_signif(data = dplyr::filter(cephalosporin, Treatment == "Control"),
mapping=aes(y = Abundance, x = Time_tx),
comparisons = list(c("Day -1","Week 1"),c("Week 5", "Week 9")),
test = "wilcox.test", step_increase = 0.04, y_position = c(19,20),
color = "#374e55", tip_length = 0.002,
test.args=list(alternative = "less", var.equal = FALSE, paired=T), textsize = 3.5) +
geom_signif(data = dplyr::filter(cephalosporin, Treatment == "Antibiotic", sample_ID %notin% "MA029.5"),
mapping=aes(y = Abundance, x = Time_tx),
comparisons = list(c("Day -1","Week 1"), c("Week 5", "Week 9")),
test = "wilcox.test", step_increase = 0.04, y_position = c(17,18),
color = "#df8f44", tip_length = 0.002,
test.args=list(alternative = "less", var.equal = FALSE, paired=T), textsize = 3.5)
cephalosporin_boxplot
Cephalosporin inactivation abundance table
cephalosporin_inactivation <- metadata %>%
dplyr::filter(`Resistance Mechanism` == "antibiotic inactivation") %>%
dplyr::select(sample_ID,`ARO Term`,normalized_abundance,Time_tx,Treatment) %>%
dplyr::rename(Allele = `ARO Term`) %>%
dplyr::filter(Allele %in% ceph$Allele) %>%
dplyr::group_by(sample_ID, Time_tx,Treatment) %>%
dplyr::summarise(Abundance = sum(normalized_abundance)) %>%
dplyr::mutate(Type="Beta-lactamases")
cephalosporin_inactivation$Treatment <- factor(cephalosporin_inactivation$Treatment, levels= c("Control","Antibiotic"))
Cephalosporin inactivation (Beta-lactamases) abundance plot
#Calculating p-values between treatments by time point
stat.test <- cephalosporin_inactivation %>%
filter(sample_ID %notin% c("MA028.7")) %>%
group_by(Time_tx) %>%
wilcox_test(Abundance ~ Treatment, alternative = "less", paired = T) %>%
adjust_pvalue(method = "bonferroni") %>%
add_significance("p.adj") %>%
add_xy_position(x = "Time_tx", dodge = 0.8)
stat.test
library(ggpubr)
library(ggsci)
betalactamases_boxplot = cephalosporin_inactivation %>%
ggboxplot(x = "Time_tx", y = "Abundance",
color = "Treatment", palette = "jama", fill = "Treatment", alpha = 0.5,
add = c("jitter"), notch = F, outlier.shape = NA,
facet.by = "Type",
xlab = F, ylab = "Abundance score") +
scale_fill_jama(alpha = 0.5) +
ylim(2,24)+
stat_pvalue_manual(stat.test, label = "p", tip.length = 0) +
geom_signif(aes(y = Abundance, x = Time_tx),
comparisons = list(c("Day -1","Week 1"), c("Week 5", "Week 9")),
test = "wilcox.test",
test.args=list(alternative = "less", var.equal = FALSE, paired=F),
color = "black", tip_length = 0.02, y_position = c(18,22)) +
geom_signif(data = dplyr::filter(cephalosporin_inactivation, Treatment == "Control"),
mapping=aes(y = Abundance, x = Time_tx),
comparisons = list(c("Day -1","Week 1"), c("Week 5", "Week 9")),
test = "wilcox.test", step_increase = 0.04, y_position = c(16.5,19.5),
color = "#374e55", tip_length = 0.002,
test.args=list(alternative = "less", var.equal = FALSE, paired=T),textsize = 3.5) +
geom_signif(data = dplyr::filter(cephalosporin_inactivation, Treatment == "Antibiotic", sample_ID %notin% "MA029.5"),
mapping=aes(y = Abundance, x = Time_tx),
comparisons = list(c("Day -1","Week 1"), c("Week 5", "Week 9")),
test = "wilcox.test", step_increase = 0.04, y_position = c(15,18),
color = "#df8f44", tip_length = 0.002,
test.args=list(alternative = "less", var.equal = FALSE, paired=T), textsize = 3.5)
betalactamases_boxplot
Table with the most abundant beta-lacatamases
beta_alleles <- metadata %>%
dplyr::rename(Allele = `ARO Term`) %>%
dplyr::filter(Allele %in% c("CfxA6","CfxA2","ACI-1")) %>%
dplyr::select(sample_ID,Allele,normalized_abundance,Time_tx,Treatment) %>%
dplyr::filter(Allele %in% ceph$Allele) %>%
dplyr::group_by(sample_ID, Time_tx,Treatment,Allele) %>%
dplyr::summarise(Abundance = sum(normalized_abundance))
beta_alleles$Treatment <- factor(beta_alleles$Treatment, levels= c("Control","Antibiotic"))
Plot of the most abundant betalactamases
#Calculating p-values between treatments by time point
stat.test <- beta_alleles %>%
group_by(Time_tx,Allele) %>%
wilcox_test(Abundance ~ Treatment, alternative = "less") %>%
adjust_pvalue(method = "bonferroni") %>%
add_significance("p") %>%
add_xy_position(x = "Time_tx", dodge = 0.8)
stat.test
alleles_boxplot = beta_alleles %>%
ggboxplot(x = "Time_tx", y = "Abundance",
color = "Treatment", palette = "jama", fill = "Treatment", alpha = 0.5,
add = c("jitter"), notch = F, outlier.shape = NA,
facet.by = "Allele", nrow = 1,ncol = 5, scales = "free",
xlab = F, ylab = "Abundance score") +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
alleles_boxplot
beta_alleles_taxa_alleles <- metadata %>%
dplyr::rename(Allele = `ARO Term`) %>%
dplyr::select(Allele,`Resistomes & Variants: Observed Pathogen(s)`) %>%
dplyr::filter(Allele %in% esbls$Allele)
head(beta_alleles_taxa_alleles)
Saving taxa associated with beta-lactamases
setwd("/Users/karlavasco/Library/CloudStorage/OneDrive-MichiganStateUniversity/Manning_lab/Mastitis_project/metagenome_KV/tables/")
write.csv(beta_alleles_taxa_alleles, "./resistome/beta_alleles_taxa_alleles.csv",row.names = T)
Aggregate phyloseq objects by time points
physeq_day1 <- subset_samples(physeq, Time_tx%in%"Day -1")
physeq_week1 <- subset_samples(physeq, Time_tx%in%c("Week 1"))
physeq_week5 <- subset_samples(physeq, Time_tx%in%c("Week 5"))
physeq_week9 <- subset_samples(physeq, Time_tx%in%c("Week 9"))
Day 1 (other time points were done in the same way)
physeq <- physeq_day1
LEfSE
lefse <- microbial::ldamarker(physeq, group="Treatment", pvalue = 0.05, normalize = T,method = "log2")
lefse_sigtab <- lefse[which(lefse$p.value<0.05), ]
lefse_sigtab <- lefse_sigtab %>% select(rank, direction,tax, LDAscore, p.value, p.adj) %>% arrange(tax)
ANCOMBC
library(nloptr)
library(ANCOMBC)
#ANCOMBC analysis comparison between treatment groups
out = ancombc(data = physeq, formula = "Treatment",
p_adj_method = "holm", lib_cut = 0,
group = "Treatment", struc_zero = F, neg_lb = FALSE,
tol = 1e-05, max_iter = 100, conserve = TRUE,
alpha = 0.05, global = TRUE)
#ANCOMBC results as a list
res = out$res
#ANCOMBC results as a table
ancom_results = res %>% as.data.frame()
# Filtering only significant results
ancom_signif <- ancom_results %>% dplyr::filter(p_val.TreatmentControl <= 0.05)
Saving Lefse and ANCOM-BC results as an Excel file
#setwd("/Users/karlavasco/Library/CloudStorage/OneDrive-MichiganStateUniversity/Manning_lab/Mastitis_project/metagenome_KV/tables/differential/microbiome_metaphlan")
#library(xlsx)
#write.xlsx(as.data.frame(lefse_sigtab), file="microbiome_day1_treatment_diff.xlsx", sheetName="lefse", row.names=FALSE)
#write.xlsx(ancom_signif, file="microbiome_day1_treatment_diff.xlsx", sheetName="ancombc", append=TRUE, row.names=FALSE)
MaAsLin2
library(Maaslin2)
setwd("/Users/karlavasco/Library/CloudStorage/OneDrive-MichiganStateUniversity/Manning_lab/Mastitis_project/metagenome_KV/tables/differential/resistome_rgi_ge")
# Formating abundance table for MaAsLin2
table <- merge(physeq@tax_table,data.frame(otu_table(physeq)), by = 0) %>% remove_rownames %>% column_to_rownames(var="Row.names")
fit_data = Maaslin2(
input_data = merge(physeq@tax_table,data.frame(otu_table(physeq)), by = 0) %>%
remove_rownames %>% column_to_rownames(var="Row.names") %>% select(6:45),
input_metadata = metadata_matrix,
output = 'maaslin2_resistome_day1_notemp',
fixed_effects = c('Treatment'),
# random_effects = c("temperature_Celsius"),
min_prevalence = 0.01,
min_abundance = 0.0,
standardize = T
)
Importing results manually curated and compared between pipelines
library(readxl)
setwd("/Users/karlavasco/Library/CloudStorage/OneDrive-MichiganStateUniversity/Manning_lab/Mastitis_project/metagenome_KV/tables/resistome/differential")
diff_alleles <- read_excel("resistome_rgi_ge_diff_summary.xlsx", sheet = "alleles")
Calculating the mean of each alelle per time-point
mean_alleles <- metadata %>% select(`ARO Term`,Time_tx,normalized_abundance) %>%
rename(Allele = `ARO Term`) %>%
group_by(Allele,Time_tx) %>% summarise(mean_allele = mean(normalized_abundance))
Making table with differential abunances including a column with mean alleles
differential <- metadata %>%
rename(Allele = `ARO Term`) %>%
merge(diff_alleles, by = "Allele") %>%
merge(mean_alleles, by = c("Allele","Time_tx"))
differential$Treatment <- factor(differential$Treatment, levels= c("Control","Antibiotic"))
Differential plot week 1
diff_week1 <- differential %>%
filter(`Week 1` >= 1,Time_tx=="Week 1") %>%
select(Allele,normalized_abundance,Treatment,mean_allele) %>%
mutate(fc = normalized_abundance/mean_allele) %>%
select(Allele,fc,Treatment) %>% mutate(time = "Week 1")
diff_week1$fc <- as.numeric(as.character(diff_week1$fc))
order_allele <- diff_week1 %>% filter(Treatment =="Antibiotic") %>%
select(Allele,fc) %>%
group_by(Allele) %>%
dplyr::summarise(mean_fc = mean(fc)) %>%
arrange(mean_fc)
week1_plot <- diff_week1 %>%
mutate(fc_flip = case_when(Treatment == "Control"
~ fc*-1,
Treatment == "Antibiotic"
~ fc*1)) %>%
ggbarplot(x="Allele", y = "fc", add="mean_se",color = "Treatment", fill = "Treatment",palette = "jama",orientation = "horiz", alpha = 0.5,legend = "right", ylab = "Mean fold change", xlab = "",facet.by = "time", order = order_allele$Allele) +
geom_hline(yintercept = c(1), linetype="dotted",
color = "black", size=1)
week1_plot
Differential plot week 5
diff_week5 <- differential %>%
filter(`Week 5` >= 1,Time_tx=="Week 5") %>%
select(Allele,normalized_abundance,Treatment,mean_allele) %>%
mutate(fc = normalized_abundance/mean_allele) %>%
select(Allele,fc,Treatment) %>% mutate(time = "Week 5")
diff_week5$fc <- as.numeric(as.character(diff_week5$fc))
order_allele <- diff_week5 %>%
filter(Treatment =="Antibiotic") %>%
select(Allele,fc) %>%
group_by(Allele) %>%
dplyr::summarise(mean_fc = mean(fc)) %>%
rbind(data_frame(
Allele=c("vanR gene in vanA cluster","mef(B)"),
mean_fc=c(0,0))) %>%
arrange(mean_fc)
week5_plot <- diff_week5 %>%
ggbarplot(x="Allele", y = "fc", add="mean_se",color = "Treatment", fill = "Treatment",palette = "jama",orientation = "horiz", alpha = 0.5,legend = "right", ylab = "Mean fold change", xlab = "",facet.by = "time", order = order_allele$Allele) +
geom_hline(yintercept = 1, linetype="dotted",
color = "black", size=1)
#+ ylim(0,2)
week5_plot
Differential plot week 9
diff_week9 <- differential %>%
filter(`Week 5` >= 1,Time_tx=="Week 9") %>%
select(Allele,normalized_abundance,Treatment,mean_allele) %>%
mutate(fc = normalized_abundance/mean_allele) %>%
select(Allele,fc,Treatment) %>% mutate(time = "Week 9")
diff_week9$fc <- as.numeric(as.character(diff_week9$fc))
order_allele <- diff_week9 %>% filter(Treatment =="Antibiotic") %>%
select(Allele,fc) %>%
group_by(Allele) %>%
dplyr::summarise(mean_fc = mean(fc)) %>%
rbind(data_frame(Allele="OXA-465",mean_fc=0)) %>%
arrange(mean_fc)
week9_plot <- diff_week9 %>%
ggbarplot(x="Allele", y = "fc", add="mean_se",color = "Treatment", fill = "Treatment",palette = "jama",orientation = "horiz", alpha = 0.5,legend = "right", ylab = "Mean fold change", xlab = "",facet.by = "time", order = order_allele$Allele) +
geom_hline(yintercept = 1, linetype="dotted",
color = "black", size=1)
#+ ylim(0,2)
week9_plot
Merging differential plots
diff_plots <- ggarrange(week1_plot,week5_plot,week9_plot,common.legend = T,ncol=3,nrow = 1, widths = c(1,1.5,1.5))
diff_plots
Final figure of cephalosporin abundance and differentially abundant alleles
beta_boxplots <- ggarrange(cephalosporin_boxplot,betalactamases_boxplot, common.legend = T, labels = c("A","B"))
beta_figure <- ggarrange(beta_boxplots, diff_plots, nrow = 2, heights = c(1,1.2), widths = c(1,1.2), labels = c("","C"))
beta_figure
Saving final figure as a PDF for further modification of allele names
pdf("/Users/karlavasco/Library/CloudStorage/OneDrive-MichiganStateUniversity/Manning_lab/Mastitis_project/metagenome_KV/figures/Fig6_cephalosporin_diff.pdf", onefile=FALSE, height = 9, width = 12)
beta_figure
dev.off()
## quartz_off_screen
## 2